;---------------------------------
; find site on which model level
;---------------------------------
function location_level(slon,slat,shgt,itim,idebug)
begin

  if(idebug.eq.1) then 
  print("") 
  print("function location_level ") 
  print("-----------------------------------------------") 
  print("lon: "+slon) 
  print("lat: "+slat) 
  print("hgt: "+shgt) 
  print("time: "+itim) 
  print("") 
  end if

  filenama = "$KAILIBDATA_ROOT/geop_T63L31.nc" 
 
  fila = addfile(filenama,"r") 

  geop = fila->var156
  lat  = fila->lat
  lon  = fila->lon

  nlat = dimsizes(lat) 
  nlon = dimsizes(lon) 
  nlev = dimsizes(geop(0,:,0,0)) 

  if(idebug.eq.1) then 
  print("nlat: "+nlat) 
  print("nlon: "+nlon) 
  print("nlev: "+nlev) 
  print("") 
  end if

; interpolate geopth to site location 
;-------------------------------------------------------------------------------------

  vtt = linint2_points (lon,lat(::-1),geop(:,:,::-1,:), True, slon, slat, 0) 

  vht = vtt(:,:,0) 

  vt1 = vht 
  vt2 = vht 

  vt1(:,0) = vht(:,0) 
  vt2(:,nlev-1) = 0.  

  if(idebug.eq.1) then 
  print("") 
  printVarSummary(vht) 
  print("") 
  end if 

  do ilev = 0, nlev - 2 
     vt2(:,ilev) = vht(:,ilev) / 2.  + vht(:,ilev+1) / 2.  
  end do 

  do ilev = 1, nlev - 1 
     vt1(:,ilev) = vht(:,ilev) /2. + vht(:,ilev-1) /2.
  end do 


  klev = nlev - 1  


  do itime = itim-1, itim-1
  do ilev  = 0, nlev - 1 
     if(shgt.gt.vt2(itime,ilev).and.shgt.lt.vt1(itime,ilev)) then 
        klev = ilev 
     end if 
  end do 
  end do 

  if(idebug.eq.1) then 
  print("") 
  print("level: "+klev) 
  print("") 
  print("function location_level done! ") 
  print("-----------------------------------------------") 
  print("") 
  end if 

return(klev)
end


;---------------------------------
; find site on which model level
;---------------------------------
function location_level_GFDL(slon,slat,shgt,itim,idebug)
begin

  if(idebug.eq.1) then 
  print("") 
  print("function location_level ") 
  print("-----------------------------------------------") 
  print("lon: "+slon) 
  print("lat: "+slat) 
  print("hgt: "+shgt) 
  print("time: "+itim) 
  print("") 
  end if

  filenama = "$KAILIBDATA_ROOT/geop_GFDL.nc" 
 
  fila = addfile(filenama,"r") 

  geop = fila->var156
  lat  = fila->lat
  lon  = fila->lon

  nlat = dimsizes(lat) 
  nlon = dimsizes(lon) 
  nlev = dimsizes(geop(0,:,0,0)) 

  if(idebug.eq.1) then 
  print("nlat: "+nlat) 
  print("nlon: "+nlon) 
  print("nlev: "+nlev) 
  print("") 
  end if

; interpolate geopth to site location 
;-------------------------------------------------------------------------------------

  vtt = linint2_points (lon,lat(::-1),geop(:,:,::-1,:), True, slon, slat, 0) 

  vht = vtt(:,:,0) 

  vt1 = vht 
  vt2 = vht 

  vt1(:,0) = vht(:,0) 
  vt2(:,nlev-1) = 0.  

  if(idebug.eq.1) then 
  print("") 
  printVarSummary(vht) 
  print("") 
  end if 

  do ilev = 0, nlev - 2 
     vt2(:,ilev) = vht(:,ilev) / 2.  + vht(:,ilev+1) / 2.  
  end do 

  do ilev = 1, nlev - 1 
     vt1(:,ilev) = vht(:,ilev) /2. + vht(:,ilev-1) /2.
  end do 


  klev = nlev - 1  


  do itime = itim-1, itim-1
  do ilev  = 0, nlev - 1 
     if(shgt.gt.vt2(itime,ilev).and.shgt.lt.vt1(itime,ilev)) then 
        klev = ilev 
     end if 
  end do 
  end do 

  if(idebug.eq.1) then 
  print("") 
  print("level: "+klev) 
  print("") 
  print("function location_level done! ") 
  print("-----------------------------------------------") 
  print("") 
  end if 

return(klev)
end
